home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / ellint.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  15.0 KB  |  491 lines
  1. /* specfunc/ellint.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author: G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_precision.h>
  26. #include <gsl/gsl_sf_ellint.h>
  27.  
  28. #include "error.h"
  29.  
  30. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  31.  
  32. inline
  33. static double locMAX3(double x, double y, double z)
  34. {
  35.   double xy = GSL_MAX(x, y);
  36.   return GSL_MAX(xy, z);
  37. }
  38.  
  39. inline
  40. static double locMAX4(double x, double y, double z, double w)
  41. {
  42.   double xy  = GSL_MAX(x,  y);
  43.   double xyz = GSL_MAX(xy, z);
  44.   return GSL_MAX(xyz, w);
  45. }
  46.  
  47.  
  48. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  49.  
  50.  
  51. /* based on Carlson's algorithms:
  52.    [B. C. Carlson Numer. Math. 33, 1 (1979)]
  53.    
  54.    see also:
  55.    [B.C. Carlson, Special Functions of Applied Mathematics (1977)]
  56.  */
  57.  
  58. /* According to Carlson's algorithm, the errtol parameter
  59.    typically effects the relative error in the following way:
  60.  
  61.    relative error < 16 errtol^6 / (1 - 2 errtol)
  62.  
  63.      errtol     precision
  64.      ------     ----------
  65.      0.001       1.0e-17
  66.      0.003       2.0e-14 
  67.      0.01        2.0e-11
  68.      0.03        2.0e-8
  69.      0.1         2.0e-5
  70. */
  71.  
  72.  
  73. int
  74. gsl_sf_ellint_RC_e(double x, double y, gsl_mode_t mode, gsl_sf_result * result)
  75. {
  76.   const double lolim = 5.0 * GSL_DBL_MIN;
  77.   const double uplim = 0.2 * GSL_DBL_MAX;
  78.   const gsl_prec_t goal = GSL_MODE_PREC(mode);
  79.   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  80.   const double prec   = gsl_prec_eps[goal];
  81.  
  82.   if(x < 0.0 || y < 0.0 || x + y < lolim) {
  83.     DOMAIN_ERROR(result);
  84.   }
  85.   else if(GSL_MAX(x, y) < uplim) { 
  86.     const double c1 = 1.0 / 7.0;
  87.     const double c2 = 9.0 / 22.0;
  88.     double xn = x;
  89.     double yn = y;
  90.     double mu, sn, lamda, s;
  91.     while(1) {
  92.       mu = (xn + yn + yn) / 3.0;
  93.       sn = (yn + mu) / mu - 2.0;
  94.       if (fabs(sn) < errtol) break;
  95.       lamda = 2.0 * sqrt(xn) * sqrt(yn) + yn;
  96.       xn = (xn + lamda) * 0.25;
  97.       yn = (yn + lamda) * 0.25;
  98.     }
  99.     s = sn * sn * (0.3 + sn * (c1 + sn * (0.375 + sn * c2)));
  100.     result->val = (1.0 + s) / sqrt(mu);
  101.     result->err = prec * fabs(result->val);
  102.     return GSL_SUCCESS;
  103.   }
  104.   else {
  105.     DOMAIN_ERROR(result);
  106.   }
  107. }
  108.  
  109.  
  110. int
  111. gsl_sf_ellint_RD_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
  112. {
  113.   const gsl_prec_t goal = GSL_MODE_PREC(mode);
  114.   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  115.   const double prec   = gsl_prec_eps[goal];
  116.   const double lolim = 2.0/pow(GSL_DBL_MAX, 2.0/3.0);
  117.   const double uplim = pow(0.1*errtol/GSL_DBL_MIN, 2.0/3.0);
  118.  
  119.   if(GSL_MIN(x,y) < 0.0 || GSL_MIN(x+y,z) < lolim) {
  120.     DOMAIN_ERROR(result);
  121.   }
  122.   else if(locMAX3(x,y,z) < uplim) {
  123.     const double c1 = 3.0 / 14.0;
  124.     const double c2 = 1.0 /  6.0;
  125.     const double c3 = 9.0 / 22.0;
  126.     const double c4 = 3.0 / 26.0;
  127.     double xn = x;
  128.     double yn = y;
  129.     double zn = z;
  130.     double sigma  = 0.0;
  131.     double power4 = 1.0;
  132.     double ea, eb, ec, ed, ef, s1, s2;
  133.     double mu, xndev, yndev, zndev;
  134.     while(1) {
  135.       double xnroot, ynroot, znroot, lamda;
  136.       double epslon;
  137.       mu = (xn + yn + 3.0 * zn) * 0.2;
  138.       xndev = (mu - xn) / mu;
  139.       yndev = (mu - yn) / mu;
  140.       zndev = (mu - zn) / mu;
  141.       epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));
  142.       if (epslon < errtol) break;
  143.       xnroot = sqrt(xn);
  144.       ynroot = sqrt(yn);
  145.       znroot = sqrt(zn);
  146.       lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
  147.       sigma  += power4 / (znroot * (zn + lamda));
  148.       power4 *= 0.25;
  149.       xn = (xn + lamda) * 0.25;
  150.       yn = (yn + lamda) * 0.25;
  151.       zn = (zn + lamda) * 0.25;
  152.     }
  153.     ea = xndev * yndev;
  154.     eb = zndev * zndev;
  155.     ec = ea - eb;
  156.     ed = ea - 6.0 * eb;
  157.     ef = ed + ec + ec;
  158.     s1 = ed * (- c1 + 0.25 * c3 * ed - 1.5 * c4 * zndev * ef);
  159.     s2 = zndev * (c2 * ef + zndev * (- c3 * ec + zndev * c4 * ea));
  160.     result->val = 3.0 * sigma + power4 * (1.0 + s1 + s2) / (mu * sqrt(mu));
  161.     result->err = prec * fabs(result->val);
  162.     return GSL_SUCCESS;
  163.   }
  164.   else {
  165.     DOMAIN_ERROR(result);
  166.   }
  167. }
  168.  
  169.  
  170. int
  171. gsl_sf_ellint_RF_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
  172. {
  173.   const double lolim = 5.0 * GSL_DBL_MIN;
  174.   const double uplim = 0.2 * GSL_DBL_MAX;
  175.   const gsl_prec_t goal = GSL_MODE_PREC(mode);
  176.   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  177.   const double prec   = gsl_prec_eps[goal];
  178.  
  179.   if(x < 0.0 || y < 0.0 || z < 0.0) {
  180.     DOMAIN_ERROR(result);
  181.   }
  182.   else if(x+y < lolim || x+z < lolim || y+z < lolim) {
  183.     DOMAIN_ERROR(result);
  184.   }
  185.   else if(locMAX3(x,y,z) < uplim) { 
  186.     const double c1 = 1.0 / 24.0;
  187.     const double c2 = 3.0 / 44.0;
  188.     const double c3 = 1.0 / 14.0;
  189.     double xn = x;
  190.     double yn = y;
  191.     double zn = z;
  192.     double mu, xndev, yndev, zndev, e2, e3, s;
  193.     while(1) {
  194.       double epslon, lamda;
  195.       double xnroot, ynroot, znroot;
  196.       mu = (xn + yn + zn) / 3.0;
  197.       xndev = 2.0 - (mu + xn) / mu;
  198.       yndev = 2.0 - (mu + yn) / mu;
  199.       zndev = 2.0 - (mu + zn) / mu;
  200.       epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));
  201.       if (epslon < errtol) break;
  202.       xnroot = sqrt(xn);
  203.       ynroot = sqrt(yn);
  204.       znroot = sqrt(zn);
  205.       lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
  206.       xn = (xn + lamda) * 0.25;
  207.       yn = (yn + lamda) * 0.25;
  208.       zn = (zn + lamda) * 0.25;
  209.     }
  210.     e2 = xndev * yndev - zndev * zndev;
  211.     e3 = xndev * yndev * zndev;
  212.     s = 1.0 + (c1 * e2 - 0.1 - c2 * e3) * e2 + c3 * e3;
  213.     result->val = s / sqrt(mu);
  214.     result->err = prec * fabs(result->val);
  215.     return GSL_SUCCESS;
  216.   }
  217.   else {
  218.     DOMAIN_ERROR(result);
  219.   }
  220. }
  221.  
  222.  
  223. int
  224. gsl_sf_ellint_RJ_e(double x, double y, double z, double p, gsl_mode_t mode, gsl_sf_result * result)
  225. {
  226.   const gsl_prec_t goal = GSL_MODE_PREC(mode);
  227.   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  228.   const double prec   = gsl_prec_eps[goal];
  229.   const double lolim =       pow(5.0 * GSL_DBL_MIN, 1.0/3.0);
  230.   const double uplim = 0.3 * pow(0.2 * GSL_DBL_MAX, 1.0/3.0);
  231.  
  232.   if(x < 0.0 || y < 0.0 || z < 0.0) {
  233.     DOMAIN_ERROR(result);
  234.   }
  235.   else if(x + y < lolim || x + z < lolim || y + z < lolim || p < lolim) {
  236.     DOMAIN_ERROR(result);
  237.   }
  238.   else if(locMAX4(x,y,z,p) < uplim) {
  239.     const double c1 = 3.0 / 14.0;
  240.     const double c2 = 1.0 /  3.0;
  241.     const double c3 = 3.0 / 22.0;
  242.     const double c4 = 3.0 / 26.0;
  243.     double xn = x;
  244.     double yn = y;
  245.     double zn = z;
  246.     double pn = p;
  247.     double sigma = 0.0;
  248.     double power4 = 1.0;
  249.     double mu, xndev, yndev, zndev, pndev;
  250.     double ea, eb, ec, e2, e3, s1, s2, s3;
  251.     while(1) {
  252.       double xnroot, ynroot, znroot;
  253.       double lamda, alfa, beta;
  254.       double epslon;
  255.       gsl_sf_result rcresult;
  256.       int rcstatus;
  257.       mu = (xn + yn + zn + pn + pn) * 0.2;
  258.       xndev = (mu - xn) / mu;
  259.       yndev = (mu - yn) / mu;
  260.       zndev = (mu - zn) / mu;
  261.       pndev = (mu - pn) / mu;
  262.       epslon = locMAX4(fabs(xndev), fabs(yndev), fabs(zndev), fabs(pndev));
  263.       if(epslon < errtol) break;
  264.       xnroot = sqrt(xn);
  265.       ynroot = sqrt(yn);
  266.       znroot = sqrt(zn);
  267.       lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
  268.       alfa = pn * (xnroot + ynroot + znroot) + xnroot * ynroot * znroot;
  269.       alfa = alfa * alfa;
  270.       beta = pn * (pn + lamda) * (pn + lamda);
  271.       rcstatus = gsl_sf_ellint_RC_e(alfa, beta, mode, &rcresult);
  272.       if(rcstatus != GSL_SUCCESS) {
  273.         result->val = 0.0;
  274.     result->err = 0.0;
  275.         return rcstatus;
  276.       }
  277.       sigma  += power4 * rcresult.val;
  278.       power4 *= 0.25;
  279.       xn = (xn + lamda) * 0.25;
  280.       yn = (yn + lamda) * 0.25;
  281.       zn = (zn + lamda) * 0.25;
  282.       pn = (pn + lamda) * 0.25;
  283.     }
  284.     
  285.     ea = xndev * (yndev + zndev) + yndev * zndev;
  286.     eb = xndev * yndev * zndev;
  287.     ec = pndev * pndev;
  288.     e2 = ea - 3.0 * ec;
  289.     e3 = eb + 2.0 * pndev * (ea - ec);
  290.     s1 = 1.0 + e2 * (- c1 + 0.75 * c3 * e2 - 1.5 * c4 * e3);
  291.     s2 = eb * (0.5 * c2 + pndev * (- c3 - c3 + pndev * c4));
  292.     s3 = pndev * ea * (c2 - pndev * c3) - c2 * pndev * ec;
  293.     result->val = 3.0 * sigma + power4 * (s1 + s2 + s3) / (mu * sqrt(mu));
  294.     result->err = prec * fabs(result->val);
  295.     return GSL_SUCCESS;
  296.   }
  297.   else {
  298.     DOMAIN_ERROR(result);
  299.   }
  300. }
  301.  
  302.  
  303. /* [Carlson, Numer. Math. 33 (1979) 1, (4.1)] */
  304. int
  305. gsl_sf_ellint_F_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result)
  306. {
  307.   double sin_phi  = sin(phi);
  308.   double sin2_phi = sin_phi*sin_phi;
  309.   double x = 1.0 - sin2_phi;
  310.   double y = 1.0 - k*k*sin2_phi;
  311.   gsl_sf_result rf;
  312.   int status = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
  313.   result->val = sin_phi * rf.val;
  314.   result->err = GSL_DBL_EPSILON * fabs(result->val) + fabs(sin_phi*rf.err);
  315.   return status;
  316. }
  317.  
  318.  
  319. /* [Carlson, Numer. Math. 33 (1979) 1, (4.2)] */
  320. int
  321. gsl_sf_ellint_E_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result)
  322. {
  323.   const double sin_phi  = sin(phi);
  324.   const double sin2_phi = sin_phi  * sin_phi;
  325.   const double x = 1.0 - sin2_phi;
  326.   const double y = 1.0 - k*k*sin2_phi;
  327.   if(x < GSL_DBL_EPSILON) {
  328.     return gsl_sf_ellint_Ecomp_e(k, mode, result);
  329.   }
  330.   else {
  331.     gsl_sf_result rf;
  332.     gsl_sf_result rd;
  333.     const double sin3_phi = sin2_phi * sin_phi;
  334.     const int rfstatus = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
  335.     const int rdstatus = gsl_sf_ellint_RD_e(x, y, 1.0, mode, &rd);
  336.     result->val  = sin_phi * rf.val - k*k/3.0 * sin3_phi * rd.val;
  337.     result->err  = GSL_DBL_EPSILON * fabs(sin_phi * rf.val);
  338.     result->err += fabs(sin_phi*rf.err);
  339.     result->err += k*k/3.0 * GSL_DBL_EPSILON * fabs(sin3_phi * rd.val);
  340.     result->err += k*k/3.0 * fabs(sin3_phi*rd.err);
  341.     return GSL_ERROR_SELECT_2(rfstatus, rdstatus);
  342.   }
  343. }
  344.  
  345.  
  346. /* [Carlson, Numer. Math. 33 (1979) 1, (4.3)] */
  347. int
  348. gsl_sf_ellint_P_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result)
  349. {
  350.   const double sin_phi  = sin(phi);
  351.   const double sin2_phi = sin_phi  * sin_phi;
  352.   const double sin3_phi = sin2_phi * sin_phi;
  353.   const double x = 1.0 - sin2_phi;
  354.   const double y = 1.0 - k*k*sin2_phi;
  355.   gsl_sf_result rf;
  356.   gsl_sf_result rj;
  357.   const int rfstatus = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
  358.   const int rjstatus = gsl_sf_ellint_RJ_e(x, y, 1.0, 1.0 + n*sin2_phi, mode, &rj);
  359.   result->val  = sin_phi * rf.val - n/3.0*sin3_phi * rj.val;
  360.   result->err  = GSL_DBL_EPSILON * fabs(sin_phi * rf.val);
  361.   result->err += n/3.0 * fabs(sin3_phi*rj.err);
  362.   return GSL_ERROR_SELECT_2(rfstatus, rjstatus);
  363. }
  364.  
  365.  
  366. /* [Carlson, Numer. Math. 33 (1979) 1, (4.4)] */
  367. int
  368. gsl_sf_ellint_D_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result)
  369. {
  370.   const double sin_phi  = sin(phi);
  371.   const double sin2_phi = sin_phi  * sin_phi;
  372.   const double sin3_phi = sin2_phi * sin_phi;
  373.   const double x = 1.0 - sin2_phi;
  374.   const double y = 1.0 - k*k*sin2_phi;
  375.   gsl_sf_result rd;
  376.   const int status = gsl_sf_ellint_RD_e(x, y, 1.0, mode, &rd);
  377.   result->val = sin3_phi/3.0 * rd.val;
  378.   result->err = GSL_DBL_EPSILON * fabs(result->val) + fabs(sin3_phi/3.0 * rd.err);
  379.   return status;
  380. }
  381.  
  382.  
  383. /* [Carlson, Numer. Math. 33 (1979) 1, (4.5)] */
  384. int
  385. gsl_sf_ellint_Kcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
  386. {
  387.   if(k*k >= 1.0) {
  388.     DOMAIN_ERROR(result);
  389.   }
  390.   else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {
  391.     /* [Abramowitz+Stegun, 17.3.33] */
  392.     const double y = 1.0 - k*k;
  393.     const double a[] = { 1.38629436112, 0.09666344259, 0.03590092383 };
  394.     const double b[] = { 0.5, 0.12498593597, 0.06880248576 };
  395.     const double ta = a[0] + y*(a[1] + y*a[2]);
  396.     const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));
  397.     result->val = ta + tb;
  398.     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  399.     return GSL_SUCCESS;
  400.   }
  401.   else {
  402.     return gsl_sf_ellint_RF_e(0.0, 1.0 - k*k, 1.0, mode, result);
  403.   }
  404. }
  405.  
  406.  
  407. /* [Carlson, Numer. Math. 33 (1979) 1, (4.6)] */
  408. int
  409. gsl_sf_ellint_Ecomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
  410. {
  411.   if(k*k >= 1.0) {
  412.     DOMAIN_ERROR(result);
  413.   }
  414.   else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {
  415.     /* [Abramowitz+Stegun, 17.3.36] */
  416.     const double y = 1.0 - k*k;
  417.     const double a[] = { 0.44325141463, 0.06260601220, 0.04757383546 };
  418.     const double b[] = { 0.24998368310, 0.09200180037, 0.04069697526 };
  419.     const double ta = 1.0 + y*(a[0] + y*(a[1] + a[2]*y));
  420.     const double tb = -y*log(y) * (b[0] + y*(b[1] + b[2]*y));
  421.     result->val = ta + tb;
  422.     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  423.     return GSL_SUCCESS;
  424.   }
  425.   else {
  426.     gsl_sf_result rf;
  427.     gsl_sf_result rd;
  428.     const double y = 1.0 - k*k;
  429.     const int rfstatus = gsl_sf_ellint_RF_e(0.0, y, 1.0, mode, &rf);
  430.     const int rdstatus = gsl_sf_ellint_RD_e(0.0, y, 1.0, mode, &rd);
  431.     result->val = rf.val - k*k/3.0 * rd.val;
  432.     result->err = rf.err + k*k/3.0 * rd.err;
  433.     return GSL_ERROR_SELECT_2(rfstatus, rdstatus);
  434.   }
  435. }
  436.  
  437.  
  438. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  439.  
  440. #include "eval.h"
  441.  
  442. double gsl_sf_ellint_Kcomp(double k, gsl_mode_t mode)
  443. {
  444.   EVAL_RESULT(gsl_sf_ellint_Kcomp_e(k, mode, &result));
  445. }
  446.  
  447. double gsl_sf_ellint_Ecomp(double k, gsl_mode_t mode)
  448. {
  449.   EVAL_RESULT(gsl_sf_ellint_Ecomp_e(k, mode, &result));
  450. }
  451.  
  452. double gsl_sf_ellint_F(double phi, double k, gsl_mode_t mode)
  453. {
  454.   EVAL_RESULT(gsl_sf_ellint_F_e(phi, k, mode, &result));
  455. }
  456.  
  457. double gsl_sf_ellint_E(double phi, double k, gsl_mode_t mode)
  458. {
  459.   EVAL_RESULT(gsl_sf_ellint_E_e(phi, k, mode, &result));
  460. }
  461.  
  462. double gsl_sf_ellint_P(double phi, double k, double n, gsl_mode_t mode)
  463. {
  464.   EVAL_RESULT(gsl_sf_ellint_P_e(phi, k, n, mode, &result));
  465. }
  466.  
  467. double gsl_sf_ellint_D(double phi, double k, double n, gsl_mode_t mode)
  468. {
  469.   EVAL_RESULT(gsl_sf_ellint_D_e(phi, k, n, mode, &result));
  470. }
  471.  
  472. double gsl_sf_ellint_RC(double x, double y, gsl_mode_t mode)
  473. {
  474.   EVAL_RESULT(gsl_sf_ellint_RC_e(x, y, mode, &result));
  475. }
  476.  
  477. double gsl_sf_ellint_RD(double x, double y, double z, gsl_mode_t mode)
  478. {
  479.   EVAL_RESULT(gsl_sf_ellint_RD_e(x, y, z, mode, &result));
  480. }
  481.  
  482. double gsl_sf_ellint_RF(double x, double y, double z, gsl_mode_t mode)
  483. {
  484.   EVAL_RESULT(gsl_sf_ellint_RF_e(x, y, z, mode, &result));
  485. }
  486.  
  487. double gsl_sf_ellint_RJ(double x, double y, double z, double p, gsl_mode_t mode)
  488. {
  489.   EVAL_RESULT(gsl_sf_ellint_RJ_e(x, y, z, p, mode, &result));
  490. }
  491.